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Abstract 



The Earth's comparatively massive moon, formed via a giant impact on the proto-Earth, has played an important role 
in the development of life on our planet, both in the history and strength of the ocean tides and in stabilizing the chaotic 
I spin of our planet. Here we show that massive moons orbiting terrestrial planets are not rare. A large set of simulations 
by Morishima et al. (2010), where Earth- like planets in the habitable zone form, provides the raw simulation data for 



P\| our study. We use limits on the collision parameters that may guarantee the formation of a circumplanetary disk after 
£>^a protoplanet collision that could form a satellite and study the collision history and the long term evolution of the 
& satellites qualitatively. In addition, we estimate and quantify the uncertainties in each step of our study. We find that 
giant impacts with the required energy and orbital parameters for producing a binary planetary system do occur with 
more than 1 in 12 terrestrial planets hosting a massive moon, with a low-end estimate of 1 in 45 and a high-end estimate 
of 1 in 4. 

Keywords: Moon, Terrestrial planets, Planetary formation, Satellites, formation 
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1. Introduction 

The evolution and survival of life on a terrestrial planet 
requires several conditions. A planet orbiting the central 
star in its habitable zone provides the temperature suit- 
able for the existence of liquid water on the surface of 
the planet. In addition, a stable climate on timescales of 
more than a billion years may be essential to guarantee a 
suitable environment for life, particularly land-based life. 
Global climate is mostly influenced by the distribution of 
solar insolati on (|Milankovitch| |1941| |Berger et al.[ |1984 



Berger 1989 Atobe and Ida[|2006[ ). The annual- averaged 
insolation on the surface at a given latitude is, beside the 
distance to the star, strongly related to the tilt of the ro- 
tation axis of the planet relative to the normal of its orbit 
around the star, the obliquity. If the obliquity is close to 
0° , the poles become very cold due to negligible insolation 
and the direction of the heat flow is poleward. With in- 
creasing obliquity, the poles get more and more insolation 
during half of a year while the equatorial region becomes 
colder twice a year. If the obliquity is larger than 57°, 
the poles get more annual insolation than the equator and 
the heat flow changes. Therefore, the equatorial region 



can even be covered by seasonal ice (Ward and Brown- 
lee 2000). Thus, the obliquity has a strong influence on a 
planet's climate. The long-term evolution of the Earth's 
obliquity and the obliquity of the other terrestrial planets 
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in the solar system, or planets in general, is controlled by 
spin-orbit resonances and the tidal dissipation due to the 
host star and satellites of the planet. Thus, the evolution 
of the planetary obliquity is unique for each planet. 

Earth's obliquity fluctuates currently ± 1.3° around 23.3 
with a period of ~ 41, 000 years ( Laskar and Robutel . 1993 



Laskar 1996 ) . The existence of a massive (or close) satel- 



lite results in a higher precession frequency which avoids 
a spin-orbit resonance. Without the Moon, the obliquity 
of the Earth would suffer very large chaotic variations. 
The other terrestrial planets in the solar system have no 
massive satellites. Venus has a retrograde spin direction, 
whereas a possibly initial more prograde spin may have 
been influenced strongly by spin-orbit resonances and tidal 
effects (IGoldreich and Pealel 119701 ILaskarl 119961). Mars' 



obliquity oscillates ±10° degree around 25° with a period 



of several 100,000 years (Ward 1974 Ward and Rudy 



1991). Mercury on the other hand is so close to the sun 
that its rotation period is in an exact 3 : 2 resonance with 
its orbital period. Mercury's spin axis is aligned with its 
orbit normal. 

On larger timescales, the variation of the obliquity can 
be even more dramatic. It has been shown that the tilt 
of Mars' rotation axis ranges from 0° to 60° in less than 
50 million years and 0° to 85° in the case of the obliquity 
of an Earth without the Moon ( jLaskar and Robutel| [l993; 



Laskar 1996) 



The main purpose of this report is to explore the giant 
impact history of the planets in order to calculate the prob- 
ability of having a giant Moon-like satellite companion, 
based on simulations done by Morishima et al. (2010). A 
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giant impact between a planetary embryo called Theia, the 
Greek titan that gave birth to the Moon goddess Selene, 



first named by Halliday (2000), and the proto-Earth is 



the accepted model for the origin of our Mo on (jHartmann 



by the last major impactor and can not be accumulated 
via the ordered accretion of small planetesimals. Giant 
impacts with a certain impact angle and velocity gener- 
ate a disk of ejected material around the target which is 
a preliminary step in the formation of a satellite. Usu- 
ally, the simulations assume perfect accretion in a colli- 
sion. Tables of the collision outcome can help to improve 
the simulations or to estimate the errors in the planetary 
lar system (Touboul et al. 2007). After its formation, the spin, (Kokubo and Genda 20101. Until recently simula- 



and Davis |1975| [Cameron and Ward} |1976[ |Cameron and 



Benz 1991 ), an event which took place within about 100 Myr 
after the formation of calcium aluminum-rich inclusions in 
chondritic meteoroids, the oldest dated material in the so- 



Moon was much closer and the Earth was rotating more 
rapidly. The large initial tidal forces created high tidal 
waves several times per day, possibly promoting the cyclic 



replication of early bio-molecules ( |Lathe 2004) and pro- 
foundly affecting the early evolution of life. Tidal energy 
dissipation has caused the Moon to slowly drift into its cur- 
rent position, but its exact orbital evolution is still part of 



an on-going debate ( Varga et al. 2006 Lathe 2006). Cal- 



culating the probability of life in the Universe ( Ward and 



Brownlee 2000 ) as well as the search for life around nearby 
planets may take into account the likelihood of having a 
massive companion satellite. 

This report is structured as follows: In section [2j we 
give a brief review on the evolution of simulating terres- 
trial planet formation with N-body codes during the last 
decades and present the method we used. In section [3j we 
study the different parameters of a protoplanet collision to 
identify potential satellite forming events. In section |1J we 
summarize the different uncertainties from the simulations 
and our analysis that may affect the final results. Finally, 
we give a conclusion, we present our results and compare 
them with previous works in section [5] 

2. Simulating terrestrial planet formation 

There are good observational data on extra-solar gas gi- 
ant planets, but whilst statistics on extra-solar rocky plan- 
ets will be gathered in the coming years, for constraints 
on the formation of the terrestrial planets we rely on our 
own solar system. The established scenario for the forma- 
tion of the Earth and other rocky planets is that most of 
their masses were built up through the gravitational colli- 



sions and interactions of smaller bodies ( 


Chamberlin 


1905 


Safranov 


1969 


Lissauer 


1993 


1. 


Wetherill and Stewart 



( 1989 ) observed the phase of run-away growth. This phase 
is characterized by the rapid growth of the largest bodies. 
While their mass increases, their gravitational cross sec- 
tion increases due to gravitational focusing. When a body 
reaches a certain mass, the velocities of close planetes- 
imals are enhanced, the gravitational focusing decreases 
and so does the accretion efficiency. This is called the oli- 



garchic growth phase, first described by (Kokubo and Ida 



1998). During this phase, the smaller embryos will grow 
faster than the larger ones. At the end, several bodies 
of comparable size are embedded in a planetesimal disk. 
These protoplanets merge via giant impacts to form the 
final planets. Dones and Tremaine ( 1993 ) showed that 



tions were limited in the number of planetesimal bodies 
that could be self-consistently followed for time spans of 
up to billions of years, but recent algorithmic improve- 



ments by Duncan et al. ( 1998 ) in his SyMBA code and by 
Chambers ( 1998[ ) in his Mercury code have allowed them 
to follow over long time spans a relatively large number of 
bodies (0(1000) ) with high precision, particularly during 
close encounters and mergers between the bodies, where 



individual orbits must be carefully integrated (Chambers 


and Wetherill 


1998 |Agnor et al. 


1999 


Raymond et al. 


2004 Kokubo et al. 


2006 




Raymond et al. ( 


2009 


) have 



also recently conducted a series of simulations where they 
varied the initial conditions for the gas giant planets and 
also track the accretion of volatile-rich bodies from the 
outer asteroid belt, leading either to "dust bowl" terres- 
trial planets or "water worlds" and everything between 
these extremes. All prior simulation methods with full in- 
teraction among all particles have however been limited in 
number of particles since their force calculations scale as 
0{N 2 ). 

We have developed a new parallel gravity code that 
can follow the collisional growth of planetesimals and the 
subsequent long-term evolution and stability of the result- 
ing planetary system. The simulation code is based on 
an O(N) fast multipole method to calculate the mutual 
gravitational interactions, while at the same time follow- 
ing nearby particles with a highly accurate mixed vari- 
able symplectic integrator, which is similar to the SyMBA 
(Duncan et al. 19981 algorithm. Since this is completely 



integrated into the parallel code PKDGRAV2 ( |Stadel[|2001| ), 
a large speed-up from parallel computation can also be 
achieved. We detect collisions self-consistently and also 
model all possible effects of gas in a laminar disk: aero- 
dynamic gas drag, disk-planet interaction including Type- 
I migration, and the global disk potential which causes 
inward migration of secular resonances with gas dissipa- 
tion. In contrast to previously mentioned studies, this 
code allows us to self-consistently integrate through the 
last two phases of planet formation with the same numer- 
ical method while using a large number of particles. 

Using this new simulation code we have carried out 64 
simulations which explore sensitivity to the initial condi- 
tions, including the timescale for the dissipation of the so- 
lar nebula, the initial mass and radial distribution of plan- 



etesimals and the orbits of Jupiter and Saturn (Morishima 



most of a terrestrial planet's prograde spin is imparted 



et al. 2010). All simulations start with 2000 equal-mass 
particles placed between 0.5 and 4AU. The initial mass 
of the planetesimal disk nid is 5 or 10 m®. The surface 
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density E of this disk and of the initial gas disk depends 
on the radius through E oc r~~ p , where p is 1 or 2. The 
gas disk dissipates exponentially in time and uniformly in 
space with a gas dissipation time scale r gas = 1,2,3 or 
5 Myr. After the disappearance of the gas disk (more pre- 
cisely after time r gas from the beginning of the simulation) , 
Jupiter and Saturn are introduced on their orbits, e.g. cir- 
cular orbits or the current orbits with higher eccentricities. 

Figure [l] shows two merger trees: a 1.8 to® planet 
formed in the simulation with (r gas ,p,md)=(l Myr, 1,5 to®) 
and gas giants on the present orbits and a 1.1 m® planet 
formed in the simulation with (T gas ,p,TOd) = (l Myr, 2, 10 to®) 
and gas giants on circular orbits. The red branches are 
satellite forming impactors both with a mass 0.3 to® and 
represent two events of the final sample in figure [7j These 
merger trees with their different morphologies reveal the 
variety of collision sequences in terrestrial planet forma- 
tion. They show that a large set of impact histories is 
generated by these simulations despite the relatively nar- 
row parameter space for the initial conditions. 



3. Satellite formation 

During the last phase of terrestrial planet formation, 
the giant impact phase, satellites form. Collisions between 
planetary embryos deposit a large amount of energy into 
the colliding bo dies and large p arts of them heat up to sev- 
eral 10 3 K, e.g. Canup (2004). Depending on the impact 

hot molten 



angle and velocity and the involved masses 
material from the target and impactor can be ejected into 
an circumplanetary orbit. This forms a disk of ejecta, the 
disk material is in a partially vapor or partially molten 
state, around the target planet. The proto-satellite disk 
cools and solidifies. Solid debris form and subsequently 
agglomerate into a satellite (|Ohtsuki| |1993| | Canup and 



Esposito 1996 Kokubo et all 2000) 



The giant impact which resulted in the Earth-Moon 



system is a very particular event ( Cameron and Bcnz 



1991 Canup 2004). The collision parameter space that 



describes a giant impact can by parametrized by 7 = 
TOj/TOtot; the ratio of impactor mass rrii to total mass in 
the collision m to t, by v = v- lnip /v csc , the impact velocity 
in units of the escape velocity v csc — \J 2Gm to t J(Xi + ?"*)> 
where r t and are the radii of target and impactor. Fur- 
thermore, it is described by the scaled impact parameter b, 
where b — indicates a head-on collision and b = 1 a graz- 
ing encounter, and the total angular momentum L. Recent 



numerical results (Canup 2008) obtained with smoothed 



particle hydrodynamic (SPH) simulations for the Moon- 
forming impact parameters require: 7 ~ 0.11, v ~ 1.1, 
b r~j 0.7 and L ~ 1.1 Lem, where Lem is the angular mo- 
mentum of the present Earth-Moon system. These sim- 
ulations also include the effect of the initial spins of the 
colliding bodies, but the explored parameter space is re- 
stricted to being close to the Moon-forming values given 
above. 





Figure 1: Two merger trees. They illustrate the accretion from the 
initial planetcsimals to the last major impactors that merge with 
the planet. Every 'knee' is a collision of two particles and the length 
between two collisions is given by the logarithm of the time between 
impacts. The thickness of the lines indicates the mass of the particle 
(linear scale). The red branch is the identified satellite forming im- 
pact in the planet's accretion history. Top: a 1.1 m® planet formed 
in the simulation with (T ga , B ,p,m^)=(l Myr,2,lOrri0) an d gas giants 
on circular orbits and a 0.3 impactor. In this case, the moon 
forming impact is not the last collision event but it is followed by 
some major impacts. Right: a 0.7 mj planet formed in the simula- 
tion with (T gas , p,m t j) = (3 Myr, 1,5 mj) and gas giants on the present 
orbits and a 0.2 impactor. It is easy to see that the moon forming 
impact is the last major impact on the planet. Although this planet 
is smaller than the upper one, it is composed out of a similar number 
of particles. In the case of the more massive disk (mj = 10), the 
initial planctesimals are more massive since their number is constant. 
Therefore, fewer particles are needed to form a planet of comparable 
mass than in the case of mj = 5. 
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If one does not focus on a strongly constrained system 
like the Earth-Moon system but just on terrestrial planets 
of arbitrary mass with satellites that tend to stabilize their 
spin axis, the parameter space is broadened. It becomes 
difficult to draw strict limits on the parameters because 
collision simulations for a wider range of impacts were not 
available for our study. Hence, based on published Moon- 
forming SPH simulations by Canup (2004, 2008), we use a 
semi-analytic expression to constrain the mass of a circum- 
planetary disk that can form a satellite. In addition, we 
include tidal evolution and study the ability of the satellite 
to stabilize the spin axis of the planet. 

3. 1 . Satellite mass and collision parameters 

We do not know the exact outcome of a protoplanet 
collision, but certainly the satellite mass is related to the 
collision parameters of the giant impact. Hence, we can 
draw a connection from these parameters to the mass of 
the final satellite. Based on the studies of the Earth's 
Moon formation, we can start with a simple scaling re- 
lation: a Mars-size impactor gives birth to a Moon-size 
satellite. Their mass ratio is mMats/^Moon ~ 10- Thus, 
to very first approximation, we can assume that an im- 
pactor mass is usually 10 times larger than the final mass 
of the satellite. Of course, this shows that only a small 
amount of material ends up in a satellite, but this state- 
ment is only valid for a certain combination of mass ratio, 
impact parameter and impact speed and usually gives an 
upper limit on the satellite mass. 

In order to get a better estimation of the satellite mass, 



we use the method obtained in the appendix of Canup 



(2008). 



There, an expression is derived that describes the mass 
of the material that enters the orbit around the target after 
a giant impact: 



?71disk 

m-tot 



/,f, pass 

m-tot 



(1) 



where the prefactor C 7 ~ 2.8(0.1/7) ' has been deter- 
mined empirically from the SPH data. m pass /rn t ot is mass 
of the impactor that avoids direct collision with the tar- 
get. It depends mainly on the impact parameter b and on 
the mass ratio 7 and can be computed by studying the 
geometry of the collision. The total impactor volume that 
collides with the target is: 



A{cj>) d<t>, 



(2) 



with 



A{ct>) = r t 2 t (<£) + r?fli(0) - Dr t sin ^ sin 0^), (3) 

where r t and r t are impactor and target radius, D = b(r{ + 
r t ) gives the distance between the centers of the bodies and 



i(4) 



cos 



D 2 



9 • 2 
ft sin 1 



,~ 1 



2_Dr,; sin < 



(4) 



D 2 



r? sin 2 ( 



2Dr t 



(5) 



Assuming a differentiated impactor with r corc ~ 0.5ri 
and repeating the above integration for this radius, the 
colliding volume of the impactor mantle is V mant \ c = Vt ~ 
V core . If we assume that the core is iron and the mantle 
dunite, the core density p C ore has roughly twice the density 
of the mantle /O man tic- The mass of the impactor that hits 
the target is m hit = PcoreKorc + /OmantieKiantic- Therefore, 
the mass that passes the target is m pass = — mhit- 

We use the full expression derived by |Canup| ( |2008[ ), 
equation ([TJ , to estimate the disk mass resulting from the 
giant impacts in our simulation. Equation |lj is correct to 
within a factor 2, if v < 1.4 and 0.4 < b < 0.7 or if v < 1.1 
and 0.4 < b < 0.8. Figure [2] illustrates the amount of 
material that is transported into orbit for the parameter 
range 0.4 < b < 0.8 (see figure B2 in |Canup| ( |2008[ ) for 
more details). It shows that a small impact parameter b 
reduces the material ejected into orbit significantly. The 
same holds for a reduction of 7, because those collisions 
are more grazing. Based of simple arguments, we use the 
limits above to identify the moon forming collision in the 
(v, 6)-plane. Details on how this assumption affects the 
result are given in sectio n [4| 



Ida et al. (1997) and Kokubo et al. (2000) studied the 



formation of a moon in a circumplanetary disk through N- 
body simulations. They found that the final satellite mass 
scales linearly with the specific angular momentum of the 
disk. The fraction of the disk material that is finally in- 
corporated into the satellite ranges form 10 to 55%. Thus, 
we assume that not more than half of the disk material is 
accumulated into a single satellite. The angular momen- 
tum of the disk is unknown and we can not use the more 
exact relationships. 

3.2. Spin- orbit resonance 

Spin-orbit resonance occurs when the spin precession 
frequency of a planet is close to one of the planet's or- 
bital precession frequencies. It causes large variation in the 



obliquity (Laskar 19961, the angle between this spin axis 



and the normal of the planet's orbital plane. An obliquity 
stabilizing satellite increases the spin precession frequency 
to a non-resonant (spin-orbit) regime. To ensure this, one 



can set a rough limit on the system parameters (Atobe et 

aT~ 



2004[ ) through 



m.. 



> 



(6) 



where m s is the mass of the satellite, m* the mass of the 
central star and a s the semi-major axis of the satellite's 
orbit and a p the semi-major axis of the planet. 

If the left term of the inequality is much larger than 
the right one, the spin precession frequency of the planet 
should be high enough to ensure that it is over the up- 
per limit of the orbital precession frequency so that spin- 
orbit resonance does not occur. Although this inequality 
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Figure 2: The disk mass resulting after a giant impact in units of the 
total mass of the colliding system relative to the impact parameter 
b based on equation |TJ. The different solid lines belong to different 
mass ratios 7. From bottom to top: 7 = 0.05,0.1,0.15,0.2,0.3 and 
0.5. These mass ratios result in different disk masses. The equation 
is valid up to a factor 2 in between the two dash lines (0.4 < b < 0.8) 
for small velocities. The disk mass is an upper bound on the satellite 



is very simplified, we try to estimate the minimum mass 
of a satellite such that it is able to stabilize the obliquity 
of its planet. 

The exact semi-major axis of a satellite after forma- 
tion is unknown but the Roche limit is the lower bound of 
its semi major axis. The Roche limit cir of the planet is 



(Murray and Dermott 1999) 



an 



3m r 



(7) 



where r s is the radius of the satellite and m s its mass and 
the mass of the planet is given by m p . Ohtsuki (1993) 



and Canup and Esposito (1996) provided detailed ana- 



lytic treatments of the accretion process of satellites in an 
impact-generated disk. Based on those studies, |Kokubo| 
et al. (20001 have shown that the true value of the ra- 



dius of satellite accretion will not diverge much from the 
Roche radius in the case of the Earth-Moon system, a 
typical satellite orbit semi-major axis in their simulations 
was a ~ 1.3 <2r. We used this approximation to estimate a 
lower bound on the satellite-planet mass ratio. We rewrite 
the Roche limit as 



OR 



( 



(8) 



where we used = (-^ p s )~ 1 m s and the fact that p p = p s 
in our model. We insert this in equation (|6| instead of a s 
and get the condition 



9m + 



(9) 



Inserting a density of 2 g cm 2 , the density of the bodies in 
the Morishima simulations, and planet semi-major axis of 



1 AU, we get a mass ratio of ~ 10 -5 , which is smaller than 
the minimum ratio of the smallest and largest particles in 
our simulations. This is a lower limit on the stabilizing 
satellite mass but the tidal evolution of the planet-satellite 
system can alter this limit dramatically. 

The orbital precession frequencies of a planet depend 
on the neighbouring or massive planets in its system. In 
the case of the Earth, Venus, Jupiter and Saturn cause 
the most important effects. To keep the spin precession 
frequency high enough, a spin period below 12 h would 



have the same effect as the present day Moon ( Laskar and 
Robutel, 19931). Hence, even without a massive satellite 



obliquity stabilization is possible as long as the planet is 
spinning fast enough. However, the moon forming impact 
provides often a significant amount of angular momentum. 
A more sophisticated analysis including the precession fre- 
quencies of the all planets involved or formed in the simu- 
lations and a better treatment of the collisions to provide 
better estimates of the planetary spins would clearly be an 
improvement but this is out of the scope of this work. 

3.3. Tidal evolution 

After its formation, even a small satellite is stabiliz- 
ing the planet's obliquity. Its fate is mainly controlled 
by the spin of the planet, the orientation of the spin axis 
of the planet relative to its orbit around the central star 
and relative to the orbital plane of the satellite and by 
possible spin-orbit resonance. Which satellites will con- 
tinue to stabilize the obliquity as they recede from the 



planet? Orbital evolution is a complicated issue (Atobe 



and Ida 2006) and there is still an ongoing debate even 



in the case of the Earth-Moon system, as mentioned. To 
classify the different orbital evolutions, it is helpful to in- 
troduce the synchronous radius, at which a circular or- 
bital period equals the rotation period of the planet. In 
the prograde satellite outwards of the synchronous 

radius will recede from the planet as angular momentum 
is tidally transferred from the planet to the satellite and 
the spin frequency of the planet decreases. In this case 
the synchronous radius will grow till it equals the satel- 
lite orbit. Angular momentum is transferred faster if the 
mass of the satellite is large, since the tidal response in the 
planet due to the satellite is greater. Large mass satellites 
will quickly reach this final co-rotation radius, where their 
recession stops. Even though their orbital radius becomes 
larger, these moons are massive enough to satisfy inequal- 
ity (|6| and avoid spin-orbit resonances. Small satellites 
will recede very slowly compared to their heavy broth- 
ers and eventually fulfill the condition ^ within the host 
star's main sequence life time. In contrast, low mass satel- 
lites that form in situ far outside the Roche limit, should 
this be possible, will probably not stabilize the spin axis. 
Intermediate mass satellites with m s ~ rayioon may re- 
cede fast enough so that they tend to lose their obliquity 
stabilizing effect during a main sequence life time. The 
Earth-Moon system shows that even in this intermediate 
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mass regime, long term stability can occur, since the Moon 
has stabilized the Earth's obliquity for billion of years. 

On the other hand, a satellite inside the synchronous 
orbit will start to spiral towards the planet, while its an- 
gular momentum is transferred to the planet. Soon, it is 
disrupted by tidal forces or will crash on the planet. How 
do we know if a satellite forms inside or outside the syn- 
chronous radius? The Roche limit Q depends on the mass 
and size of the bodies while the synchronous radius r sync 
is a function of the planet mass m p and depends inversely 
proportional on its rotation frequency, which is obtained 
from equating the gravitational acceleration and the cen- 
tripetal acceleration: 



15 



r syn c = i.Gni p ) 3 uj p 3 . 
Equating this formula with ^ gives: 



-V. 



(Gm p )3tVp 



(10) 



(11) 



The planet mass drops out and we get an lower limit on 
the planet angular velocity to guarantee a satellite outside 
the Roche radius: 



J P , min = ^G Pp Y- = 0.00043 s- 1 , 



(12) 



which equals a rotation period of 4 h. The final planets of 



Morishima et al. (2010) have generally a high rotational 
speed, some of them rotating above break up speed. The 
rotation period after the moon forming collision is usu- 
ally around 2-5 hours, see figure [3j A more exact parti- 
cle growth model without the assumption of perfect stick- 
ing would lower the rotation frequency by roughly 30% 
( |Kokubo and Genda| |2010 ), where we can also include 
the loss of rotational angular momentum due to satellite 
formation (full circles). On the other hand, the mean dis- 
tance of satellite formation is around 1.3 aR, and the max- 
imum rotation period of the planet for a receding satellite 
changes to ~ 6 h (dashed line) . Hence, applying this to 
our final sample, 1/4 of all moon forming collisions are 
excluded. 

Finally, for the remaining events, we assume that the 
planet spin is large enough so that the synchronous radius 
is initially smaller than the Roche radius. Satellites that 
form behind the Roche limit will start to recede from the 
planet. We can conclude that almost every satellite-planet 
system in our simulation will fulfill ([6|. 

A special outcome of the tidal evolution of a prograde 



planet-satellite system is described by Atobe and Ida ( 2006 ) 



If the initial obliquity 6 of the planet after the moon form- 
ing impact is large, meaning that the angle between planet 
spin axis and planet orbit normal is close to 90°, results 
in a very rapid evolution when compared to the previ- 
ously discussed case of a moon receding to the co-rotating 
radius and becoming tidally locked. The spin vectors of 
the protoplanets are isotropically distributed after the gi- 
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Figure 3: The mass of the planet in units of its final mass 
mtot /m t ot. final relative to the rotational period in days of the body 
after the moon-forming collision. This is the final sample but with- 
out excluding moon formation inside the synchronous radius. Empty 
circles: The decrease of the planet spin due to escaping material or 
the formation of a satellite are not included in the calculation of the 
rotation period, it is based on the perfect accretion assumption. Full 
circles: The angular momentum is assumed to be 30% smalle r due 
to a more realistic collision model (Kokubo and Gcnda 2010). The 
dotted line gives the position of the threshold for synchronous rota- 
tion for the Roche radius cir, the dashed line gives its position for 
the 1.3 ctR. We want to focus on the better estimate of the anuglar 
momentum (full circles) and the initial moon radius (dashed line). 
Hence, roughly a quarter of all moon forming collisions are excluded 
in the final sample, figure [7] 



distribution corresponds to p(9) = \ sin(t9) . In the most 
extreme scenario, a massive prograde satellite will crash 
onto the planet in a timescale of order 10,000 years after 
formation, even though it initially recedes. Hence, without 
a favorable initial obliquity is an added requirement for the 
survivability of a close or massive moon. We include this 
by discarding massive and highly oblique impacts by an 
approximated inequality, based on area B in figure 13 in 



Atobe and Ida (20061 



< 



7T m s 
0.2 m„ 



(13) 



ant impact phase (Agnor et al. 1999) and the obliquity 



More exactly, this holds best for planets with a p ~ 1 AU, 
and it has only to be taken into account if m s /m p < 0.05. 
The distribution of the a p of the simulated planets ranges 
for 0.1 to 4 AU. In the case of a smaller distance to the host 
star, angular momentum is removed faster from the planet- 
satellite system and the evolution timescales are shorter in 
general. To use this limit properly, a more general expres- 
sion for different semi-major axes has to be derived, which 
is out of the scope of this work. Moreover, it depends 
linearly on the satellite mass which is overestimated in 
general. A smaller satellite mass would reduce the num- 
ber of excluded events in general. However, when applying 
this constraint on the data, only one out of seven moon 
forming collisions are affected. Since it is over-simplified, 
we exclude it from our analysis. 

For retrograde impacts, where the impactor hits in op- 
position to the target's spin, two cases result in differing 
evolution. If the angular momentum of the collision is 
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much larger than the initial rotational angular momen- 
tum, the spin direction of the planet is reversed and any 
impact-generated disk will rotate in the same direction as 
the planet's spin. On the other hand, a retrograde collision 
with small angular momentum will not alter the spin di- 
rection of the planet significantly and it becomes possible 



lost during collisions in general, but as mentioned above, 



to be left with a retrograde circumplanetary disk(Canup 



2008). After accretion of a satellite, tidal deceleration due 



to the retrograde protoplanet spin will reduce the orbital 
radius of the bodies continuously till they merge with the 
planet. Hence, a long-lived satellite can hardly form in 
this case. 

In order to find a reasonable threshold between these 
two regimes based on the limited information we have, we 
assume that the sum of the initial angular momenta of the 
bodies L t and Li and of the collision L co \ equals the spin 
angular momentum of the planet and the satellite £ p i ane t 
and Lmoon after the collision plus the angular momentum 
of the orbiting satellite £ or bit at 1.3 Or parallel to the col- 
lision angular momentum. 

Comparing initial and final angular momenta gives: 



L C a\ + L t + Li — L 



orbit 



L 



planet 5 



where £ rbit is parallel to L co \. 

-^orbit — |iorbit|"T=i . • 

l^coll 



(14) 



(15) 



We assume that |L m oon|/|(£orbit + ^planet )| < 1 and 



Gi 



iorbit = m s a 2 s n = m s a 2 s ^ = m s ^Ga s m p , (16) 

where n is the orbital mean motion of the satellite if its 
mass is much smaller than the planet mass, n ~ */ Guip/af, 
and G is the gravitational constant. With a s — we get 



3 1 

iorbit = m s {Gm p )i { TOj 

\2tt p p 



If 



l-^col + L t \ 



< 1, 



the collision is retrograde. Hence, 

-^planet — -^col + Lt + Li — |L rbit 



L 



col 



\L, 



(17) 



(18) 



(19) 



col 



If ipianct and L ibit, which is parallel to L C oi, are retro- 
grade, 



l-^col + L 

orbit | 
V( L lol + ^orbit) 



<1, 



(20) 



the satellite will be tidally decelerated. Those cases are 
excluded from being moon-forming events. Spin angular 
momenta are in general overestimated since material is 



the simulations of Morishima et al. (2010) assume perfect 
accretion. We include this consideration by reducing the 
involved spins and the orbit angular momentum of the 



satellite by 30%, ( |Kokubo and Genda| [2010| . 

Both scenarios described above, a large initial obliquity 
or a retrograde orbiting planet, might not become impor- 
tant until subsequent impactors hit the target. Giant im- 
pacts can change the spin state of the planet in such a 
way that the satellite's fate is to crash on the planet. This 
scenario is discussed in the next section. 

3.4- Collisional history 

We exclude all collisions from being satellite-forming 
impacts whose target is not one of the final planets of a 
simulation. A satellite orbiting an impact is lost through 
the collision with the larger target. 

Multiple giant impacts occur during the formation pro- 
cess of a planet and it is useful to study the impact his- 
tory in more detail. Subsequent collisions and accretion 
events on the planet after the satellite-forming event may 
have a large effect on the final outcome of the system. 
We choose a limit of 5TOpi an ctcsimai to distinguish between 
large impacts and impacts of small particles, which are 
responsible for the ordered accretion. To stay consistent, 
the same limit is used below to exclude small impactors 
from our analysis. We divide all identified moon forming 
events into four groups: 

a) The moon forming event is the last major impact on 
the planet. Subsequent mass growth happens basi- 
cally through planetesimal accretion (see merger tree 
at the bottom in figure [T]). 

b) There are several moon forming impacts, in which 
the last impact is the last major impact on the planet. 

c) The moon forming event is not the last giant im- 
pact on the planet. The satellite can be lost due to 



a disruptive near or head-on-collision (Stewart and 



Leinhardt 2009) of an impactor and the satellite. A 



late giant impact on the planet can change the spin 
axis of the planet and the existing satellite can get 
lost due to tidal effects (see tree at the top in figure 

d) There are several moon forming events in the impact 
history of the planet, followed by additional major 
impacts. As before, the moon forming collisions can 
remove previously formed satellites. On the other 
hand, and existing moon can have an influence on 
the circumplanetary disk formed by a giant impact 
and can suppress the formation of multiple satellites 
orbiting the planet. 

The final states of the planet-satellite systems in group c 
are difficult to estimate, since such systems might change 
significantly by additional giant impacts. To a lesser ex- 
tent, this holds for d, but those systems are probably more 
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Figure 4: This bar chart shows the distribution of the moon forming 
collisions in four groups with different impact history with respect 
to the last major impact, a: the moon forming impact is the last 
major impact on the planet, b: there are multiple moon forming 
impacts, but the last one is also the last major impact, c: there 
is only one moon forming impact and it is followed by subsequent 
major impacts, d: there are multiple moon forming impacts, but the 
last of them is followed by subsequent major impacts. 



resilient to the loss of satellites by direct collisions. The 
number of events per group is shown in figure [4j Group c 
and d include more then 2/3 of all collisions. Group c is 
the most uncertain and we use it to quantify the error on 
the final sample. 

Furthermore, we exclude impactors and targets that 
have masses of the order of an initial planetesimal mass 
(5 mpianetesimai ~ 0.0025 m®) from producing satellite form- 
ing events since their masses are discretized and related to 
the resolution of the simulation. In addition, small im- 
pactors will probably not have enough energy to eject a 
significant amount of material into a stable orbit. There- 
fore, setting a lower limit on the target and impactor mass 
results in the exclusion of many collisions but not in a sig- 
nificant underestimation of the true number of satellites, 
see figure [5j 

4. Uncertainties 

Our final result depends on several assumptions, lim- 
itations and approximations. In this section we want to 
quantify them as much as possible and summarize them. 



The data of Morishima et al. (20101 we are using has 



two peculiarities worth mentioning: The focus on the Solar 
System and the small number of simulations per set of 
initial conditions. 

The simulations were made in order to reproduce the 
terrestrial planets of the Solar System. The central star 
has 1 solar mass and the two gas giants that are intro- 
duced after the gas dissipation time scale have the mass 
of Jupiter and Saturn and the same or similar orbital el- 
ements. Although the initial conditions like initial disk 
mass or gas dissipation time scale are varied in a certain 
range, the simulated systems do not represent general sys- 
tems with terrestrial planets. However, we assume that 
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Figure 5: The target mass rat versus the impactor mass m; in- 
volved in the satellite-forming collisions. We exclude events that 
include impactors and targets with masses of the order of the ini- 
tial planetesimal mass {rai, rat < 5 m p i anctes ; ma i) indicated by the 
dashed line to avoid resolution effects. The line in the case of the 
target mass is not shown since it is very close the frame. Due to this 
cut, the total number of satellites decreases significantly while the 
number of massive satellites in our analysis increases. The shown 
sample is the final set of events without applying the threshold for 
small particles. Therefore, the number of accepted events in this 
plot does not equal the number of satellites in figure [7] since this cut 
is not applied. If the cut is used, new events are accepted, which 
where neglected before because they were not the last moon forming 
impacts on the target. 



the range of impact histories is representative of the range 
that would be seen in other systems. Merger trees (figure 
[TJ reveal that those simulations cover a huge diversity of 
impact histories. But future work will need to investigate 
the full range of impact histories that could be relevant to 
the formation of terrestrial planets in other, possibly more 
exotic, extra-solar systems. 

However, the set of simulation covers a broad range 
of initial conditions. But for every set of initial condi- 
tions, only one simulation exists since they are very time 
consuming (each simulation requires about 4 months of 
a quad-core CPU). Therefore, it is difficult to separate ef- 
fects of the choice of certain initial parameters from effects 
of stochastic processes. We grouped our moon forming 
events with respect to the simulation parameter in ques- 
tion. The only parameter that reveals an effect on the final 
sample is the gas dissipation time scale T gas . Larger time 
scales lead to less moon forming collisions. This variation 
is correlated with mass and number of final planets. If 
the gas disk stays for several million years, the bodies are 
affected by the gas drag for a long time, spiral towards 
the sun and get destroyed. Therefore, there are less giant 
impacts and smaller planets. One would suppose that the 
initial mass of the gas disk should be correlated to the mass 
of the final planets and therefore to the number of giant 
impacts, but the two initial protoplanetary disk masses of 
5 m® and 10 m® show no essential difference. 

The approach we use to identify events and estimate 
the mass of the satellite is also based on various approxi- 
mations and limitations. 
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Satellite mass. The method we use to calculate the mass 
of the circumplanetary disk is valid to better than a factor 



of 2 within the parameter range we use (Canup 2008) 



Only 10-55% of the disk mass are embodied in the satellite 



(Kokubo et al. 2000). Therefore, the mass we use is just 



half of the estimated disk mass, and in the worst case, the 
satellite is ten times less massive than estimated. This 
uncertainty affects the number of massive satellites but 
not the number of satellites in general. 

In the (v,b)-plane, we use the same restrictive limits to 
constrain the collision events. Figure [6] presents all moon 
forming events and this parameter range. Hence, this area 
gives just a lower bound. We see that it covers some of 
the most populated parts, but a significant amount of the 
collisions are situated outside this area. Equation ([I]) can 
help to constrain the outcome of the collisions close outside 
of the shaded region. A collision with a small impact pa- 
rameter (b < 0.4) will bring very little material into orbit, 
whatever the impactor mass involved. In this regime, we 
have few events with high impact velocity and even with 
high velocity it might be very hard to eject a significant 
amount of material into orbit. Hence, we exclude them 
from being moon forming events. For intermediate impact 
parameter (0.4 < b < 0.8), there are high velocity events 
(v > 1.4). Above a certain velocity threshold, depending 
on b and 7, most material ejected by the impact will es- 
cape the system and the disk mass might be too small to 
form a satellite of interest. A large parameter (b > 0.8) 
describes a highly grazing collision. It is difficult to extent 
equation ([I]) for larger 6, since these collisions will proba- 
bly result in a hit-and-run events for high velocities. SPH 
simulations ( |Canup| 2004, 2008) suggest that high rela- 
tive impact velocities (v > 1.4) will increase rapidly the 
amount of material that escapes the system. Neverthe- 
less, these sets of simulations focus not on general impacts 
and the multi-dimensional collision parameter space is not 
studied well enough to describe those collisions in more 
detail. Detailed studies of particle collisions will hopefully 



be published in the near future (e.g. Kokubo and Genda 



(2010)). Based on the arguments above, the events inside 



of the shaded area form our final sample. Most of the colli- 
sions outside the area will not form a moon. The collisions 
with 0.4 < b < 0.8 and velocity slightly above v = 1.4 and 
with b > 0.8 and small velocity (v < 1.4) are events with 
unknown outcome. Including those events, the final num- 
ber of possible moon forming events is increased by not 
more than a factor of 1.5. 

Tidal evolution - Rotation period. In order to separate re- 
ceding satellites from satellites which are decelerated after 
their formation, we check if the initial semi-major axis 
of the satellite is situated outside or inside of the syn- 
chronous radius of the planet. In figure [3j the sample is 
plotted twice, once including a general correction for an- 
gular momentum loss due to realistic collisions and once 
without correction. Moreover, two thresholds for the syn- 
chronous radius are shown. We choose the threshold at 



1.3ap{ (dashed line) and the corrected rotation period (full 
circles) to be the most justified case. There are two most 
extreme cases: the threshold situated at Ior (dotted line) 
and a corrected rotation period (full circles) indicates that 
only one in five collisions lead to receding moons. On the 
other hand, the threshold situated at 1.3 csr and a rota- 
tion period directly obtained form the simulations (empty 
circles) indicates that only one in eight collisions lead to a 
non-receding moon. 

Tidal evolution - Retrograde satellites. To exclude retro- 
grade orbiting satellites, we use a simple relation between 
the angular momenta involved in a collision. Since we have 
no exact data of the angular momentum distribution af- 
ter the collision, this limit is very approximate. To get 
an estimate of the quality of this angular momentum ar- 
gument, we study the effect of this threshold on the final 
sample. First, roughly half of all moon-forming impacts 
are retrograde. But in almost every case, L rbit is much 
smaller than the initial angular momentum. Only four of 
the retrograde collisions have not enough angular momen- 
tum to provide a significant change of the spin axis and 
those planet-satellite systems remain retrograde. There- 
fore, our final estimate is not very sensitive to this angular 
momentum argument. A larger set of particle collision 
simulations could provide better insight into the angular 
momentum distribution. 

Collision history. Two issues affecting the collision history 
can change the final result. These being the mass cut we 
choose to avoid resolution effects, where we exclude small 
impactors and targets, and the uncertainty about the fi- 
nal state of the planet-satellite system because of multiple 
and subsequent major impacts. The first limit seems to be 
well justified. If we change the threshold mass or exclude 
the cut completely, we change mainly the number of very 
small, non-satellite forming, impacts since in these cases 
there is usually insufficient energy to bring a significant 
amount of material into orbit. The second issue has a sig- 
nificant effect on the final sample, however. Group c, the 
single moon forming collisions followed by large impacts, 
contains the most uncertain set of events: Neither all satel- 
lites survive the subsequent collisions nor is it likely that 
all satellites get lost. Hence, in the extreme case, almost 
half of all initially formed moons, all of Group c, could be 
lost and our final result reduced significantly. 

An overview of the above uncertainties in given in ta- 
ble 1. Additional planet formation simulations are nec- 
essary to quantify a large part of existing uncertainties. 
Simulations of protoplanet collisions exploring the multi- 
dimensional collision parameter space are desirable and 
will hopefully be published soon and might help to con- 
strain the parameters for moon formation better. A study 
of the effect of subsequent accretion of giant bodies after 
a moon forming impact might give better insight in the 
evolution of a planet-satellite system. 
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Figure 6: The parameter space in the (b, t>j mp /?; esc )-plane. The 
darker region gives the region for with equation Q holds up to a 
factor of 2. The shown sample is the final sample without the con- 
straints in the (b, i>)-plane. Similar to figure [5] the number of ac- 
cepted events will increase when applying the cut, since an earlier 
moon forming impact on the planet can become suitable. As one 
expects, there are more collisions with larger impact parameter. Ig- 
noring events with smaller b or higher v (see text), there is still a 
significant group of events with b > 0.8. Much material escapes in 
the high velocity cases and no moon will form, but the low velocities 
arc hard to exclude or calculate. 



Table 1 : A list of the different conditions that affect the final number 
of satellites. The uncertainty factor gives the range in which the 
final result varies around the most justified. I equals the final value. 
The first two factors can not be estimated, it shows that additional 
simulations would be helpful. The estimation of the satellite mass 
is separated from the rest of the list, since the uncertainty on this 
estimate does not change the number of satellites in the final sample, 
in contrast to the others. The exclusion of retrograde satellites is 
very approximate, but it has almost no effect on the final sample 
and therefore this factor is close to 1. 



5. Discussion and results 

Under these restrictive conditions we identify 88 moon 
forming events in 64 simulations, the masses of the re- 
sulting planet-satellite systems are shown in figure [7] On 
average, every simulation gives three terrestrial planets 
with different masses and orbital characteristics and we 
have roughly 180 planets in total. Hence, almost one in 
two planets has an obliquity stabilizing satellite in its or- 
bit. If we focus on Earth-Moon like systems, where we 
have a massive planet with a final mass larger than half 
of an Earth mass and a satellite larger than half a Lunar 
mass, we identify 15 moon forming collisions. Therefore, 
1 in 12 terrestrial planets is hosting a massive moon. The 
main source of uncertainties results from the modelling of 
the collision outcomes and evolution of the planet-satellite 
system as well as the small number of simulation and the 
limited range of initial conditions. We do not include the 
latter in our estimate. Hence, we expect the total number 
of Earth- Moon like systems in all our simulations to be in 
a range from 4 to 45. This results in a low-end estimate 
of 1 in 45 and a high-end estimate of 1 in 4. In addition, 
taking into account the uncertainties on the estimation of 
the satellite mass, roughly 60 of those systems are formed 
in the best case or almost no such massive satellites are 
formed if the efficiency of the satellite accretion in the cir- 
cumplanetary disk is very low. 

There are several papers, where the authors performed 
N-body simulations and searched for moon-forming colli- 
sions. 



Agnor et al. (1999) started with 22-50 planetary 
embryos in a narrow disk centered at 1 AU. They esti- 
mated around 2 potentially moon-forming collisions per 
simulation, where the total angular momentum of the en- 
counter exceeds the angular momentum of the Earth-Moon 
system. They pointed out that this number is somewhat 
sensitive to the number, spacing, and masses of the ini- 
tial embryos. O'Brien et al. (20061 performed simula- 
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tions with 25 roughly Mars-mass embryos embedded in 
a disk of 1000 non-interacting (with each other) planetes- 
imals in an annulus from 0.3 to 4.0 AU. They found that 
giant impact events which could form the Moon occur 
frequently in the simulations. These collisions include a 
roughly Earth-size target whose last large impactor has a 
mass of 0.11 — 0.14Mb and a velocity, when taken at in- 



finity, of 4km/s, as found by Canup (2004). O'Brien et al 



pointed out that their initial embryo mass is close to the 



impactor mass. Raymond et al. ( 2009 ) set up about 90 em- 



bryos with masses from 0.005 to O.IMe in a disk of more 
than 1000 planetesimals, again with non-interaction of the 
latter. Assuming again Canup's requirements (v/v esc < 
1.1, 0.67 < shi9 < 0.76, 0.11 < 7 < 0.15), only 4% of 
their late giant impacts fulfill the angle and velocity cri- 
teria. They concluded that Earth's Moon must be a cos- 
mic rarity but a much larger range of late giant collisions 
would produce satellites with different properties than the 
Moon. The initial embryo size seems to play a role in 
those results. In contrast, the simulations of Morishima et 
al. (2010) start with 2000 fully interacting planetesimals 



and since embryos form self-consistently out of the plan- 
etesimals in these simulations, problems with how to seed 
embryos are completely avoided. To constrain the simula- 



tions, Morishima et al. (2010) were interested in the tim- 
ing of the Moon- forming impact. To identify potentially 
events, they searched for a total mass of the impactor and 
the target > O.5m0, a impactor mass > O.O5to0 and a im- 
pact angular momentum > Lem- They found almost 100 
suitable impacts in the 64 simulations. Since their sample 
also includes high velocity or grazing impacts, although 
constrained through the more general angular momentum 
limit, and does not take into account collision history and 
tidal evolution, the difference to our result is not surpris- 
ing. 

Life on planets without a massive stabilizing moon 
would face sudden and drastic changes in climate, posing 
a survival challenge that has not existed for life on Earth. 
Our simulations show that Earth-like planets are common 
in the habitable zone, but planets with massive, obliquity 
stabilizing moons do occur only in 10% of these. 
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